
####################################################
####          Replication Code for              ####
####     Bayesian versus maximum likelihood     ####
####      estimation of treatment effects       ####
####          in bivariate probit               ####  
####       instrumentalvariable models          ####    
####################################################
####                                            ####  
####                4/ 18/2017                  ####         
####     this file runs rcode to estimate       ####
####      the Stan models                       ####
####################################################



#FORMAT DATA FROM Sondheimer and Green#
####################################################################

##perry model 
perry.data = read_dta("perry.data.dta")
N = dim(perry.data)[1]
dat = list(K=1,D=2,N=N,y = as.matrix(data.frame(perry.data$diplomaonly,perry.data$antvote)), x = matrix(rep(1,N),ncol=1),instr = c(perry.data$treatment), SD =2.5)
rstan_options(auto_write = TRUE)
options(mc.cores = 4)
perry.fit = stan(file = "BivarProbit.stan", data=dat, iter=10000, chains=4,thin=5, verbose=FALSE, seed = 12345, control= list(adapt_delta = 0.95))
save(perry.fit, file = "perry.fit.rda")


### ihad model
ihad.data = read_dta("ihad.data.dta")
N = dim(ihad.data)[1]
dat = list(K=1,D=2,N=N,y = as.matrix(data.frame(ihad.data$diplomaonly,ihad.data$antvote)), x = matrix(rep(1,N),ncol=1),instr = c(ihad.data$treatment), SD = 2.5)
ihad.fit = stan(file = "BivarProbit.stan", data=dat, iter=10000, chains=4,thin=5, verbose=FALSE, seed = 12345, control= list(adapt_delta = 0.95))
save(ihad.fit, file ="ihad.fit.rda")


###star model
star.data = read_dta("star.data.dta")
N =dim(star.data)[1]
dat = list(K=1,D=2,N=N,y = as.matrix(data.frame(star.data$diplomaonly,star.data$antvote)), x = matrix(rep(1,N),ncol=1),instr = c(star.data$treatment), SD = 2.5)
star.fit = stan(file = "BivarProbit.stan", data=dat, iter=10000, chains=4,thin=5, verbose=FALSE, seed = 12345, control= list(adapt_delta = 0.95))
save(star.fit, file ="star.fit.rda")


